#####################################
#### Assign necessary libraries. ####

starttime <- Sys.time()

library(survival)
library(plyr)
library(optmatch)
library(rcbalance)
library(rcbsubset)
library(car)
library(data.table)
library(splitstackshape)
library(RCurl)
library(MASS)
library(pROC)


source("/Users/sharpe/Desktop/Young Surgeons/code/apply_penalty.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/subtabfun.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/meantab_both.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/baltab_ys.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/smahal_func.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/penalize.near.exact.v2.R")

###############################################
#### Read in the match and risk datasets. ####


ortho_risk_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/ortho_risk_patients.csv")
risk_score_patients_trad <- subset(ortho_risk_patients, era_trad==1)
risk_score_patients_mod <- subset(ortho_risk_patients, era_mod==1)

ortho_match_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/ortho_match_patients.csv")
trad_match_patients <- subset(ortho_match_patients, era_trad==1)
mod_match_patients <- subset(ortho_match_patients, era_mod==1)

NEW_PATIENTS_BEFORE <- subset(ortho_match_patients, new_surgeon_flag==1 & era_mod==1)
exp_patients_before <- subset(ortho_match_patients, new_surgeon_flag==0 & era_mod==1)

###############################################
###############################################

ortho_npgs <- c("hip_repair", "hip_replacement", "hip_replacement_revision", "knee_replacement", "knee_replacement_revision")
ortho_npgs_for_prop_score <- c("hip_repair", "hip_replacement", "hip_replacement_revision", "knee_replacement")


comorbidities <- c('silb_AbdCancer',  'silb_ANGINA', 'silb_ASTHMA', 'silb_CANCER', 'silb_CHF', 'silb_CHRNLUNG', 'silb_ChrPepticUlcer',
                   'silb_COAG', 'silb_CollagenVascular', 'silb_CongenitalCoag', 'silb_CUSHINGS', 'silb_DEMENT', 'silb_DIABETES', 'silb_GRAVES', 'silb_HTN', 'silb_HYPOTHY', 'silb_LIVER',
                   'silb_PARAPLEG', 'silb_PastARR', 'silb_PastMi', 'silb_PostPulmFibr', 'silb_RENLDYS', 'silb_renlfail',
                   'silb_SEIZUR', 'silb_SickleCell', 'silb_STROKE', 'silb_THROMB', 'silb_UnstAngina', 'silb_Valvulardis', 'silb_AIDS')
aids_position <- match("silb_AIDS", comorbidities)
#### Remove AIDS from the comorbidities, since it is pretty rare.
comorbidities_for_prop_score <- comorbidities[-aids_position]

prop_score_pat_vars_mod <- c("new_surgeon_flag","admsn_in_quarters", ortho_npgs_for_prop_score, "emergent_type")


prop_score_pat_mod_post_2010 <- glm(formula = new_surgeon_flag ~ ., data=mod_match_patients[which(mod_match_patients$oct_2010_plus==1),c(prop_score_pat_vars_mod)], family=binomial)
prop_score_pat_mod_pre_2010 <- glm(formula = new_surgeon_flag ~ ., data=mod_match_patients[which(mod_match_patients$oct_2010_plus==0),c(prop_score_pat_vars_mod)], family=binomial)

summary(prop_score_pat_mod_post_2010)
summary(prop_score_pat_mod_pre_2010)


mod_match_patients$prop_prob_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="link", link="logit"),0)

mod_match_patients$prop_prob_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="link", link="logit"),0)

#### Calculate the standard deviation for the propensity scores. ####

prop_score_sd_mod_post_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==1), c("prop_log_pat_mod_post_2010")])
prop_score_sd_mod_pre_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==0), c("prop_log_pat_mod_pre_2010")])

###############################################################################
#### Create a variable to prevent matrix/dataframe from becoming a vector. ####

mod_match_patients$for_matrix <- 1



#############################
#### Match the patients. ####


mod_match_patients <- mod_match_patients[order(-c(mod_match_patients$treatment)),]

mod_treat <- mod_match_patients$treatment


maha_vars_pat_imp <- c("admsn_in_quarters", ortho_npgs)

mod_match_patients$exact_group_2 <- paste(mod_match_patients$pair_id_surg, mod_match_patients$oct_2010_plus, mod_match_patients$Narrow_Procedure_Group, sep="_")



exact_vars_pat <- c("exact_group_2")
caliper_var_pat_mod <- c("admsn_in_quarters", "for_matrix", "prop_log_pat_mod_pre_2010", "prop_log_pat_mod_post_2010", "emergent_type")

near_exact_vars <- NA



mnf.names <- vector("list",length=1)
mnf.names[[1]] <- c("admit_year", "emergent_type")

maha_data_pat_imp_mod <- mod_match_patients[,maha_vars_pat_imp]
exact_matrix_pat_mod <- mod_match_patients[,exact_vars_pat]
calip_data_pat_mod <- mod_match_patients[,caliper_var_pat_mod]



################################################################
#### Create the important distance matrices; add penalties. ####

dist_struc_pats_imp_mod <- build.dist.struct(z=mod_treat, X = maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, calip.option="none")

avg_imp_mod <- mean(unlist(dist_struc_pats_imp_mod))

for(i in 1:length(dist_struc_pats_imp_mod)){
  for(j in 1:length(dist_struc_pats_imp_mod[[i]])){
    dist_struc_pats_imp_mod[[i]][j] <- (dist_struc_pats_imp_mod[[i]][j]/avg_imp_mod)*1500
  }
  names(dist_struc_pats_imp_mod[[i]]) <- names(dist_struc_pats_imp_mod[[i]])
}

max(unlist(dist_struc_pats_imp_mod))


####

dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_post_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_post_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_pre_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_pre_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="emergent_type")

####

dist_struc_pats_mod <- dist_struc_pats_imp_mod

####


treated_info_pats_mod <- subset(mod_match_patients, treatment==1)
control_info_pats_mod <- subset(mod_match_patients, treatment==0)

patient_match_mod <- rcbsubset(distance.structure = dist_struc_pats_mod, fb.list = mnf.names, treated.info = treated_info_pats_mod, control.info = control_info_pats_mod,
                               penalty=2, tol=10)


tem_match_pats_mod <- treated_info_pats_mod[as.numeric(rownames(patient_match_mod$matches)),]
sda_match_pats_mod <- control_info_pats_mod[(patient_match_mod$matches),]
pair_id_pats_mod <- c(seq(1:dim(tem_match_pats_mod)[1]), seq(1:dim(sda_match_pats_mod)[1]))
combined_pats_mod <- rbind(tem_match_pats_mod, sda_match_pats_mod)
combined_pats_mod$pair_id_pats <- pair_id_pats_mod

#### Add propensity scores to the before datasets.

NEW_PATIENTS_BEFORE$prop_prob_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)


exp_patients_before$prop_prob_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="link", link="logit"),0)

NEW_PATIENTS_BEFORE$treatment <- 1



###############################################
#### Save necessary pieces of information. ####

#### Save the before match datasets. ####
write.csv(x=NEW_PATIENTS_BEFORE, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_ortho_new_pats.csv")
write.csv(x=exp_patients_before, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_ortho_exp_pats.csv")



#### Save the matched data. ####
write.csv(combined_pats_mod,  "/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_ortho_matched_mod_pats.csv")

#### Save the models. ####
prop_score_mod_summary <- summary(prop_score_pat_mod_pre_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_2_propensity_model_mod_ortho_pre_2010.csv")
prop_score_mod_summary <- summary(prop_score_pat_mod_post_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_2_propensity_model_mod_ortho_post_2010.csv")



######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################

rm(list=ls())

######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################

#####################################
#### Assign necessary libraries. ####

starttime <- Sys.time()

library(survival)
library(plyr)
library(optmatch)
library(rcbalance)
library(rcbsubset)
library(car)
library(data.table)
library(splitstackshape)
library(RCurl)
library(MASS)
library(pROC)


source("/Users/sharpe/Desktop/Young Surgeons/code/apply_penalty.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/subtabfun.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/meantab_both.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/baltab_ys.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/smahal_func.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/penalize.near.exact.v2.R")

###############################################
#### Read in the match and risk datasets. ####


gensurg_risk_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/gensurg_risk_patients.csv")
risk_score_patients_trad <- subset(gensurg_risk_patients, era_trad==1)
risk_score_patients_mod <- subset(gensurg_risk_patients, era_mod==1)

gensurg_match_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/gensurg_match_patients.csv")
trad_match_patients <- subset(gensurg_match_patients, era_trad==1)
mod_match_patients <- subset(gensurg_match_patients, era_mod==1)

NEW_PATIENTS_BEFORE <- subset(gensurg_match_patients, new_surgeon_flag==1 & era_mod==1)
exp_patients_before <- subset(gensurg_match_patients, new_surgeon_flag==0 & era_mod==1)


################################
#### Define some variables. ####


gensurg_npgs <- c("cholecystectomy", "hernia_abdominal_open", "colectomy_partial__open", "appendectomy",                               
                  "lysis_of_adhesions", "mastectomy", "proctopexy",
                  "hernia_groin_open", "small_bowel_other", "pyloroplasty",
                  "proctectomy", "stomach_antireflux", "large_bowel_other",
                  "thyroidectomy_partial", "biliary_other", "small_bowel_resection",
                  "pd_access_procedure", "ulcer", "pancreatectomy",
                  "parathyroidectomy", "esophagectomy", "adrenalectomy",
                  "stomach_gastric_bypass_nonbariatric", "stomach_other", "splenectomy",
                  "liver_partial_hepatectomy", "colectomy_total__open", "stomach_partial_gastrectomy",
                  "thyroidectomy_total", "stomach_total_gastrectomy", "bariatric",
                  "biliary_common_duct", "liver_other", "thyroidectomy_substernal",
                  "esophagomyotomy", "hernia_abdominal_lap", "colectomy_partial__lap",
                  "colectomy_total__lap", "hernia_groin_lap")

gensurg_npgs_for_prop_score <- gensurg_npgs[gensurg_npgs != "colectomy_partial__open"]

comorbidities <- c('silb_AbdCancer',  'silb_ANGINA', 'silb_ASTHMA', 'silb_CANCER', 'silb_CHF', 'silb_CHRNLUNG', 'silb_ChrPepticUlcer',
                   'silb_COAG', 'silb_CollagenVascular', 'silb_CongenitalCoag', 'silb_CUSHINGS', 'silb_DEMENT', 'silb_DIABETES', 'silb_GRAVES', 'silb_HTN', 'silb_HYPOTHY', 'silb_LIVER',
                   'silb_PARAPLEG', 'silb_PastARR', 'silb_PastMi', 'silb_PostPulmFibr', 'silb_RENLDYS', 'silb_renlfail',
                   'silb_SEIZUR', 'silb_SickleCell', 'silb_STROKE', 'silb_THROMB', 'silb_UnstAngina', 'silb_Valvulardis', 'silb_AIDS')
aids_position <- match("silb_AIDS", comorbidities)
#### Remove AIDS from the comorbidities, since it is pretty rare.
comorbidities_for_prop_score <- comorbidities[-aids_position]

prop_score_pat_vars_mod <- c("new_surgeon_flag","admsn_in_quarters", gensurg_npgs_for_prop_score, "emergent_type")


prop_score_pat_mod_post_2010 <- glm(formula = new_surgeon_flag ~ ., data=mod_match_patients[which(mod_match_patients$oct_2010_plus==1),c(prop_score_pat_vars_mod)], family=binomial)
prop_score_pat_mod_pre_2010 <- glm(formula = new_surgeon_flag ~ .-stomach_total_gastrectomy-biliary_common_duct-liver_other-proctopexy, data=mod_match_patients[which(mod_match_patients$oct_2010_plus==0),c(prop_score_pat_vars_mod)], family=binomial)

summary(prop_score_pat_mod_post_2010)
summary(prop_score_pat_mod_pre_2010)

mod_match_patients$prop_prob_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="link", link="logit"),0)

mod_match_patients$prop_prob_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="link", link="logit"),0)

#### Calculate the standard deviation for the propensity scores. ####

prop_score_sd_mod_post_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==1), c("prop_log_pat_mod_post_2010")])
prop_score_sd_mod_pre_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==0), c("prop_log_pat_mod_pre_2010")])

###############################################################################
#### Create a variable to prevent matrix/dataframe from becoming a vector. ####

mod_match_patients$for_matrix <- 1

#############################
#### Match the patients. ####

mod_match_patients <- mod_match_patients[order(-c(mod_match_patients$treatment)),]

mod_treat <- mod_match_patients$treatment


maha_vars_pat_imp <- c("admsn_in_quarters", gensurg_npgs, "emergent_type")

mod_match_patients$exact_group_2 <- paste(mod_match_patients$pair_id_surg, mod_match_patients$oct_2010_plus, sep="_")



exact_vars_pat <- c("exact_group_2")
caliper_var_pat_mod <- c("admsn_in_quarters", "for_matrix", "prop_log_pat_mod_pre_2010", "prop_log_pat_mod_post_2010", "emergent_type")

near_exact_vars <- NA

more_import <- c("mastectomy","thyroidectomy_partial", "thyroidectomy_total","bariatric")
middle_import <- c("adrenalectomy", "colectomy_total__open","colectomy_total__lap", "esophagectomy", "liver_partial_hepatectomy", "pancreatectomy", "proctectomy", "stomach_total_gastrectomy")
less_import <- c(gensurg_npgs[!(gensurg_npgs %in% c(more_import, middle_import, "liver_transplant"))])


mnf.names <- vector("list",length=4)
mnf.names[[1]] <- c(more_import)
mnf.names[[2]] <- c(mnf.names[[1]], middle_import)
mnf.names[[3]] <- c(mnf.names[[2]], less_import)
mnf.names[[4]] <- c(mnf.names[[3]], "admit_year", "emergent_type")


maha_data_pat_imp_mod <- mod_match_patients[,maha_vars_pat_imp]
exact_matrix_pat_mod <- mod_match_patients[,exact_vars_pat]
calip_data_pat_mod <- mod_match_patients[,caliper_var_pat_mod]



################################################################
#### Create the important distance matrices; add penalties. ####

dist_struc_pats_imp_mod <- build.dist.struct(z=mod_treat, X = maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, calip.option="none")

avg_imp_mod <- mean(unlist(dist_struc_pats_imp_mod))

for(i in 1:length(dist_struc_pats_imp_mod)){
  for(j in 1:length(dist_struc_pats_imp_mod[[i]])){
    dist_struc_pats_imp_mod[[i]][j] <- (dist_struc_pats_imp_mod[[i]][j]/avg_imp_mod)*1500
  }
  names(dist_struc_pats_imp_mod[[i]]) <- names(dist_struc_pats_imp_mod[[i]])
}

max(unlist(dist_struc_pats_imp_mod))


####

dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_post_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_post_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_pre_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_pre_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="emergent_type")

dist_struc_pats_mod <- dist_struc_pats_imp_mod



treated_info_pats_mod <- subset(mod_match_patients, treatment==1)
control_info_pats_mod <- subset(mod_match_patients, treatment==0)


patient_match_mod <- rcbsubset(distance.structure = dist_struc_pats_mod, fb.list = mnf.names, treated.info = treated_info_pats_mod, control.info = control_info_pats_mod,
                               penalty=2, near.exact=c(gensurg_npgs), tol=10)




tem_match_pats_mod <- treated_info_pats_mod[as.numeric(rownames(patient_match_mod$matches)),]
sda_match_pats_mod <- control_info_pats_mod[(patient_match_mod$matches),]
pair_id_pats_mod <- c(seq(1:dim(tem_match_pats_mod)[1]), seq(1:dim(sda_match_pats_mod)[1]))
combined_pats_mod <- rbind(tem_match_pats_mod, sda_match_pats_mod)
combined_pats_mod$pair_id_pats <- pair_id_pats_mod

#### Add propensity scores to the before datasets.

NEW_PATIENTS_BEFORE$prop_prob_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)



exp_patients_before$prop_prob_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="link", link="logit"),0)

NEW_PATIENTS_BEFORE$treatment <- 1




###############################################
#### Save necessary pieces of information. ####

#### Save the before match datasets. ####
write.csv(x=NEW_PATIENTS_BEFORE, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_gensurg_new_pats.csv")
write.csv(x=exp_patients_before, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_gensurg_exp_pats.csv")



#### Save the matched data. ####
write.csv(combined_pats_mod,  "/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_gensurg_matched_mod_pats.csv")

#### Save the models. ####
prop_score_mod_summary <- summary(prop_score_pat_mod_pre_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_2_propensity_model_mod_gensurg_pre_2010.csv")
prop_score_mod_summary <- summary(prop_score_pat_mod_post_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_2_propensity_model_mod_gensurg_post_2010.csv")





